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We compute 0{a'^lna) corrections to the decay rates of para- and orthopositronium into two 
and three photons, respectively. For this calculation we employ the nonrelativistic QED regularized 
. . . dimensionally and we explain how in this framework the logarithms of the fine structure constant 

' can be extracted. 
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I. INTRODUCTION 



Positronium decays into two and three photons provide an interesting test of bound states Quantum Electrodynamics 
(QED). While the parapositronium (p-Ps) decay rate is well described by QED, it is known that the decay of 
orthopositronium (o-Ps) into three photons is still a controversial issue Q. However, on the theoretical side there 
I , has been a major breakthrough recently and both decay rates are now known with 0{a^) accuracy (2]j|]. The next 
level of precision, i.e. 0{a^) correction, is currently beyond the reach, although parts of the 0{a^) correction that 
are enhanced by the logarithms of the fine structure constant, can be computed. 

The knowledge of these corrections is not of significant phenomenological importance at present since they are much 
' smaller than the experimental accuracy for both p-Ps and o-Ps decays. Nevertheless, from a theoretical viewpoint, it 
. is an interesting problem because at order there are both leading 0{a'^ In^ a) and subleading 0{a'^ In a) corrections. 
00 ' The calculation of the leading 0{a^ In^ a) corrections is a fairly simple enterprise; it has been done quite some time 
ago 1^. Here we are interested in subleading logarithmic corrections. In order to compute them, we utilize the 
nonrelativistic QED in dimensional regularization and we explain how in this framework the logarithmic corrections 
can be computed using limited amount of information. 
f~| Before diving into the description of the calculation, let us summarize our results. For the 0{a^ In a) corrections 

Qh^ to para- and orthopositronium decay rates we find: 

&: Arp = |^l„a|-^ + 10I„2-2^,}r<»>, (1) 

d ^ The coefficients Ap^o describe O^a/n) corrections to the lowest order decay widths Pp^o. They are [^j6| 

Ap = — - 5, Ao = -10.286606(10). (3) 

Numerically, the 0{a^lna) corrections, Eqs.(|^,^), evaluate to 7.7919 a'^/nlna Pj?'' for p-Ps and to 
—5.517 a'^/iTlna Po°'' for o-Ps. These corrections are therefore quite comparable with the "leading" 0{a'^ log^ a) 
corrections computed in 

The paper is organized as follows. In the next Section we set up the framework of the calculation. We then continue 
with detailed discussion of how various contributions to Ps decays at 0{a^ In a) are computed. In the last Section we 
present our conclusions. 
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II. FRAMEWORK OF THE CALCULATION 



Let us first discuss the framework of our calculation. We work in nonrelativistic QED regularized dimensionally; 
d = 3 — 2e is the number of spatial dimensions and e is the regularization parameter. General features of this technique 
have been described at length in our previous paper Here we would like to discuss a new issue which was not 
considered in |^ - how logarithmic in a corrections can be extracted. 

In order to extract logarithms of the fine structure constant in a self-consistent way, we use the fact that the matrix 
element of any operator in dimensional regularization is a uniform function of the fine structure constant. This 
implies that all the dependence on the fine structure constant can be scaled out of any dimensionally regularized 
matrix element^. The scaling, however, should be done in d dimensions. 

In order to establish the scaling rules, we need to know how different quantities involved in bound state calculations 
scale with a. To do that, we rewrite the d-dimensional Schrodinger equation in "atomic units" familiar from the 
standard treatment of hydrogen atom in three-dimensional Quantum Mechanics. 

Consider the Schrodinger equation for positronium in d dimensions: 



p c{d)a 
m 7-^-2 



= E^, (4) 



where c{d) = T{d/2— V) /ir'^/'^^^ (see Q). Let us rescalep fry and r r/j and choose 7 = {ma/2) ^+^' . Then both 
the Coulomb Hamiltonian H = / m — c{d)a / r'^~'^ and its eigenvalue E scale as /m. The wave function of a bound 
state, being normalized to unity, scales as 7^^/^. Hence, when expressed in atomic units, the energy of a bound state 
and its properly normalized wave function depend on d only. The scaling rules above provide sufficient information 
to write the matrix element of any operator in atomic units and therefore scale out all the dependence on the fine 
structure constant. 

There is another important point that makes the extraction of the In a corrections possible using limited amount of 
information. In order to explain it, we remind the reader that in bound state calculations different contributions to 
the final result can be distinguished. In particular, there are so-called hard contributions. The 0{a^) hard corrections 
to the Ps decay rates are described by the three-loop Feynman diagrams for the process e+e^ 2(3)7, that have to 
be computed exactly at the threshold. Schematically, such diagrams generate the correction 



^hard - (^) 



4 + - + <3 



£2 e 



VBorn, (5) 



( ii_3 are some constants) to the annihilation kernel Vboih that is responsible for the lowest order decay rate: 

^^pl = (*I^Bor„|*) (X K (6) 

Here stands for the positronium wave function at the origin. The effective potential (|5|), in turn, generates the 
correction to the decay rate. 



^hardi 



(*|Vhard|*) CX 



h t2 



£2 e 



K (7) 



If we rewrite V^o in atomic units, Eq.(|^) generates logarithms of the fine structure constant. These logarithms are 
artificial, since we anticipate that other, soft scale, contributions also generate divergences which exactly match and 
cancel all the divergences in Vhard- This implies, that the logarithms associated with the rescaling of the wave function 
at the origin get cancelled as well. Therefore, the easiest way to avoid considering T4iaid (which is not available at 
present) is to work with relative, rather than absolute, corrections to the decay width. This automatically discourages 
Eq.(^ as the source of logarithms of the fine structure constant since in this case there is simply nothing to rescale. 

In order to illustrate how these arguments help to compute the In a corrections, let us consider the matrix element of 
a nonrelativistic operator O that delivers 0{a^) correction to the lowest order annihilation kernel VBorn- In accordance 
with the above comment we consider relative correction to the decay rate: 



^We stress that this is the feature of dimensional regularization and it is not valid in other regularization schemes. 
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The operator O is a function of coordinate and momentum operators that act on the positronium wave function. 
After rescahng of all the quantities on the right hand side of Eq.(||) as described above, we end up with the following 
equation: 



(^'ll^Bornl*) 



(9) 



where n = 1, 2 is a power of a that explicitly enters O and I, k are some integer numbers. If the matrix element is 
finite, we can safely put e = and then the relative correction to the decay width is times the a-independent 
ratio and hence no logarithms of a appear. Therefore, after the rescaling, the only place where In a can come from 
is the expansion of the factor 'y3~"+''= in powers of e; this implies that in order to generate the In a corrections, the 
nonrelativistic matrix elements should diverge and only divergent pieces of the matrix elements have to be known to 
determine logarithmic in a corrections to the decay rate. Note also, that because of the relation between 7 and a. 
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even the operators that scale as integer powers of 7 can generate In a corrections. 

Let us note that the extraction of the In a corrections in dimensional regularization can lead to some counter- 
intuitive results; for example, logarithms of the fine structure constant are generated by the operators, that in a more 
"physical" regularization schemes, such as e.g. the schemes that use either the photon mass or the momentum cut-off 
A, can only lead to the logarithms of m/A but not to In a. Therefore, it appears that individual contributions to the 
final result are scheme dependent. Nevertheless, we would like to stress that once dimensional regularization and clear 
rules for extracting In a are adopted, there is no other way as to consider all possible contributions; none of them can 
be disregarded by invoking the fact that in a different regularization scheme a particular operator cannot generate 
©(In a) correction. 

The calculation of the nonrelativistic contributions that are relevant at order 0{a^ In a) is described in the following 
Sections. Some useful integrals, that we need in the calculation, are summarized in Appendix. We give intermediate 
formulas for the relative logarithmic corrections to the positronium decay rate expressed in units of (Q;'^/7r)r(l -|- 
e)3(47^)3^ 

III. IRREDUCIBLE CONTRIBUTION 

This particular contribution arises as the average value of a local operator that comes from the Taylor expansion of 
the one-loop corrections to the annihilation kernel in spatial momenta p of electron and positron. For our purpose we 
need 0{ap'^) correction to T^orn- Since this operator is constructed by Taylor expanding in external momentum, the 
non-analytic dependence on cannot appear. Then, from the rescaling argument we know that only divergent piece 
of the Wilson coefficient of this operator is required. It turns out that this divergent piece can easily be computed 
using rather general arguments. According to the rules of nonrelativistic QED, we extract Wilson coefficients of 
various effective operators from the corresponding on-shell scattering amplitudes. The key observation is that the 
divergence in the Wilson coefficient of the 0{ap'^) operator causes the divergence in the on-shell annihilation process 
e"'"e~ n"f, n = 2,3, and that it is in fact a "true" infrared divergence that should be compensated by the real 
emission of an additional soft photon in the same process. 

In order to avoid confusion we stress that the mechanism of cancelling this divergence by real radiation does not 
apply to the bound state because of C-parity conservation. There is no contradiction, however. The real infrared 
divergences in the bound state calculations are absent because electron and positron in positronium are off shell. The 
divergence appears only when we put them on mass shell in order to extract the Wilson coefficient of the relevant 
operator. The crucial observation is that considering the process in a different kinematic regime (on-shell annihilation), 
we easily find a divergent piece of the appropriate Wilson coefficient from the known amplitude of the real soft photon 
emission. 

Requiring that virtual and real corrections to the on-shell annihilation process e~^e~ — > 717 add up to a finite 
quantity, we get the following correction to the annihilation kernel Vboih: 

^irr = f^V^orn. (10 

Applying the rescaling arguments, we end up with the following correction induced by the irreducible operator ([l^): 

Ai„-^lna. (11) 
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IV. "HARD LOOP" CONTRIBUTIONS 

These contributions arise in the second order of the nonrelativistic perturbation theory. This means that the 
corresponding nonrelativistic operators O are of the form VGV, where G is the reduced Green function of the 
Coulomb Hamiltonian from Eq. and V, V are some local operators with one of the two originating from a hard 
one-loop correction. 

There are two sources of the "hard loop" contributions. The first one is the one-loop renormalization of the 
annihilation kernel: 

Vp,o = -^p,o^Born- (12) 
TT 

The corresponding 0{a^ In a) correction to the decay rate is then easily related to the 0(0^ In a) correction computed 
in We find: 

Ap,o= Qs2-2j^p,olna. (13) 
The second "hard loop" contribution corresponds to the "hard" piece in 0{ma^) effective potential. It reads fl^: 



r(l + £) 
'3to2 (47r)-^ 



ion / 32 

- + — -12 In 2 + ^+6 In 2 I - 2 In TO 



e 5 V 3 



S{r), (14) 



where S is the operator of the total spin. There is no problem with defining total spin operator here, since it multiplies 
explicitly finite quantity. 

The 0{a^) correction to the decay rate generated by the potential from Eq.(|l^) then reads: 

6uT ^ 2 mVBornGVu\-i') , (15) 

and is proportional to the Green function at the origin, G(0, 0). All necessary results for this Green function can be 
found in 0. Finally, we obtain: 

Ay = + 1 + 2 In 2 - I - 5^ (^y + In 2) + In m) In a. (16) 

Let us note that since there is an explicitly divergent term in Eq.(^), one may wonder whether or not the difference 
in the energy E of the bound state in d and three dimensions should be taken into account. We have checked that 
the cancellation of all divergent terms in the final result for the 0{a^ \na) correction to the decay rate occurs for 
arbitrary E and for this reason we use three dimensional expression for this (rescaled) quantity to present individual 
contributions as well. 

V. SEAGULL CONTRIBUTION 

The correction to the decay rate, 

5,r-2(*|VBornGV;|*), (17) 

is induced by the double seagull effective potential Vs, 

m2 J (2^)^* kk'{k + k') ' ^ ' 

where k' — q — k and Vij{k) — Sij — kikj/k'^. 

To facilitate the calculation of this potential, we introduce auxiliary integration variable fco and rewrite Vs{q) as 
follows: 



r°° dka f d'^k , , , 



(19) 
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with 



fc2 + fc2 ■ 



It is convenient to consider Fourier transform of Vs{q): 

, N f°° dkn _ , , _ , , 



where Tij{r) stands for 



To compute Tijir) we use 



{2ttY k'^ + kl kl 



(27r) 



47r 



fc2 A:2 + fc2 



d k_^,kr 47r 



(27r)'' 



2 K._,ikor), 



(20) 



(21) 



(22) 



(23) 



where K^{x) is the modified Bessel function of the second kind. 

Inserting ( psf ) to (p^, and integrating over fco in Eq.(|l]), we obtain: 



27rm2 



(47r) 



We then rescale the relative correction to the decay rate and obtain: 



5j: _ a^Y V{'^- 



*(r) 



(24) 



(25) 



As we explained previously, we need only a divergent part of the integral in Eq. (^5|) . Since in three dimensions the 
Green function G(r,0) behaves as 0{r~^) for small values of r, we may expand the wave function in series around 
r = and keep only two first terms in such an expansion: 



^ ^ = 1 - — — + 0{r^^^'^). 



4-d 



(26) 



The Green function G can be written as a sum of three pieces G = Go + Gi + Gmuiti according to the number of 
Coulomb interactions^. In three dimensions, Go ~ r"^, G\ ~ Inr and Gmuiti ~ r° as r ^ 0. Therefore, the first term 
from the expansion Eq. ( p6| ) is sufficient to extract the singularities caused by the contributions of G\ and Gmuiti ■ For 
Gmuiti we derive: 



d'^r r-3+4^Gmuiti(r-, 0) = 4^ / dr r-i+2^Gmuiti(0, 0) + 0{\) 

Jo 



Oil). 



In order to analyze Go,i contributions, it is convenient to switch to the momentum space. We obtain: 

4^7r3/2-'r(e) f d<^p Gi(p) 



d'^r r-3+4^Gi(r,0) = , , . , , o 
^ ' ' r(3/2-2e) J {2'kY P^" 



ne 



l+2c 



Oil) 



(27) 

(28) 
(29) 



The integrals in Eqs.(E3,E9|) can be expressed through the integrals listed in Appendix: 



Go is the free Green function, Gi is the single Coulomb correction to Go, and Gmuiti accounts for two and more Coulomb 
interactions. 
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-2/o(l,z). (31) 



Substituting all the relevant expressions into Eq.(p5|) and expanding in e, we arrive at the final result for the seagull 
contribution: 

Sln^a /I 41n2 5 \, 



VI. RETARDATION 



In this Section we discuss the retardation effect, caused by the exchange of a photon with a typical momentum of 
the order of the inverse Bohr radius k ^ ma between electron and positron in the bound state. The corresponding 
formula reads: 



Jretr = 2 ( * 



^BornG 



(*|VBorn|*)(^' 



Ji e 



(27r) 



2k k + H 



E 



■P - 



ikrr 



2k (k + H-E) 



■P - 

2 



ikr^ 



H.c. 



(33) 



where j'^ — p'^/m — [erk, (T]/(4m) and = pP/m + [cr'k, cr']/(4m) are electron and positron currents, respectively. 

The appearance of two matrix elements in Eq. ( p3| ) is related to the fact that we need the second order correction to 
the wave function of a bound state. In contrast to the first order, at second order of perturbation theory one should 
be careful to maintain the normalization of the wave function. This is the reason why the second term in Eq.(|3^) 
appears. The general formula for the second order correction to the wave function can be found in fl^ . 

When writing Eq. ( |33| ) , we have assumed that a magnetic photon with the momentum k is emitted by the electron 
at the point re and absorbed by the positron at the point Vp. Using the fact that H — E ^ ma^ and k ^ ma for the 
retardation effects, we expand the r.h.s. of Eq.(p3|) in powers of {H — E)/k and obtain: 



^,etr = 2 (^-ll/BornGK-etl*) - (* | ■^^Born | 



(*lK'etl*> 



where the "retardation potential" V^ct reads: 



(34) 



(35) 



and = dVret/dE. Let us illustrate how one deals with such expressions using the spin-dependent parts of the 



currents as an example. The corresponding contribution to Vrot then reads: 



T/spm _ 
''ret — 



16m^ 



d'^fc [(Tk,ai][cr'k,a'^] 
{2t:Y 1^ 



{{H-EY 



ikr 



ikr 



{H-E)^+R.c. 



(36) 



where the relative coordinate r — — Vp has been introduced. It is now easy to see that if we insert this result into 
Eq.(|^ and use the Schrodinger equation both for the wave function 4" and for the reduced Green function G, we 
obtain a finite correction to the decay rate. As we explained in the Introduction, finite contributions cannot generate 
logarithms of the fine structure constant. 

The analysis of the spin-independent contribution to Vi-ot is only slightly more cumbersome. We again pull out 
H — E both to the right and to the left by commuting it with either electron or positron current and obtain the 
following expression for the retardation potential: 



Vrct — ~- 



d'^k Vijik) 



m2 J (2t:Y fc3 



{[h-E, e^fc^p.p,} + [p., [H^p,]^^'') + H.c. 



where {H — E, e^^^piPj} denotes the anticommutator of the two operators. Using this expression for the retardation 
potential in Eq.(|3l) and applying equations of motion, we arrive at the following correction to the decay rate: 
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„d-2 



•^-t^ = (*I^Bor„G [C/y [C,p^] ,pj] I*) , (37) 



where C — ~c{d)a/r is the Coulomb potential, and Uij{r) is defined as 

Performing rescaling and computing the double commutator in Eq.(^), we get the result for the relative correction 
to the decay rate, 

...r_ 8«,--r(2-.)r(3/2-.) /,._-3.4.^(,^o)^^ (39^ 



r(o) 3m2 7r3/2-'^ J *o 

1 p.o ^ 

which is very similar to the correction induced by the seagull potential, Eq.(p5|). We can, therefore, borrow much of 
the analysis from the previous Section. We finally obtain: 

o / 2 4 In 2 \ 

A,.ct = 21nQ;- h4-41nTO Ina. (40) 

V 3e 3 / 



VII. ULTRASOFT CONTRIBUTION 



By definition, the ultrasoft contribution is due to the photons with energy and momentum of order ma^. Such soft 
photons cannot resolve the structure of the bound state and for this reason they interact directly with the positronium 
as a whole. Since the positronium is chargeless, the interaction is necessarily of the dipole nature. 

A general formula for the correction to the decay rate caused by ultrasoft contributions is: 



SusT = 2 ( * Vsorn G ^ 

\ m 

-(*|143orn|*>/* 



d'^k 47ra 'P^j{k) 2pj 
{2nY 2k E-H-k 'm 



1p^ 



d'^k ina Vi-j{k) 2pj 
{277^ 2k {E-H-kf ~^ 



After integration over directions of k and performing the rescaling we get: 



5us^ 
p(0) 

-L p.o 



l2-^7ra72-4^ nd d-l 



{2ttY d 



Gp 



dk k 



k + H 



[H,p] 













Jo 



dk fci 



{k + H - E)' 



P 



(41) 



(42) 



Both matrix elements in Eq.(^2|) must be computed in atomic units; also \^ — ^E'/^'o, and fid — 27r'^/2^^/r(3/2 — e) is 
the d-dimensional angular volume. 

The second matrix element in Eq.(^) is easy to compute^: 



P 



dk k 



l-2e 



P 



(* |p(i? - £;)-2'^p| _ 1 



2e 



{k + H-Ef 

Then, consider the first matrix element from Eg. (p^ . It is convenient to write it as a sum: 

dk 



Gp 







H -E 



[H,P] 



(43) 



(44) 



where the two terms correspond to the number of Coulomb interactions between the moments of emission and 
absorption of the ultrasoft photon^. 



^Recall that only the terms that are singular for e — > are needed. 

*If two or more Coulomb photons are exchanged, the resulting matrix element becomes finite and, in accordance with the 
argument given in Section II, it cannot generate ©(Ina) correction. 
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Let us consider Mq. In this case, integrating over k we obtain: 



G (Ho-E) ^'p [H,p] 



where Hq = /2 is the free Hamiltonian in atomic units. Using 

p [H,p] = [p, [H,p]] + [H,p] p = ATrS{r) 
we represent Mq as the sum of two terms, 

Mo = Ms + M^, 

where 



[H,p]p, 



Ms 



47r 
1 

2^ 



GiHo~E) 



-2e 



GiHo-Ey^' [H,p]p 



(45) 

(46) 
(47) 

(48) 
(49) 



We now analyze the two terms in Eq.(^7|) separately. In order to compute Ms we use the fact that only two first 
terms in the expansion of the Green function in the number of Coulomb exchanges, Go and Gi, diverge at zero spatial 
separation. Since Gniuiti(0,0) is finite, we can compute it for d — On the other hand, both Go and Gi are known 
explicitly (see e.g. |^) and hence the corresponding integrals can easily be computed. Using expressions for integrals 
summarized in Appendix, we arrive at the following expression for Ms: 



Ms 



/o(l + 2e,0) 



-/i(l + 2e,l,l)--. 



(50) 



Consider the second matrix element M^. In this case, it is sufficient to separate G = Go + (G — Gq). Since there is 
an overall factor in Eq.(^9|), the finite matrix element that contains G — Go can be calculated in three dimensions. 
Using the expression for G- Go from |], we derive: 



Aff -l(o|(G-Go)[H,p] p|*) = ^. 
The matrix element containing Go can be rewritten in momentum space as 



2e 



Go {Ho-Ey [H,p]p 



1 
2^ 



|2d 



/2 



1 



47r(p' - p)p vl/(p) 

{p'-py ~W 



(51) 



(52) 



where ^'(p) is the Fourier transform of 'l'(r). Since has an overall divergence and hence we need the integral up 
to a constant only, we can use the three dimensional expression for the wave function ^'(p): 



4:77^ r d'^p'd'^p 

~ J 



2d 



p. 



1 



ip' - p)p 
{p' - pY 



p^ 



1 



(53) 



We now rewrite the scalar product in the numerator as a linear combination of denominators and obtain our final 
expression for = + M^: 



M^ = [/i(2e, 1, 2) - /i(l + 2e, 1, 1) - /i(l + 2e, 0, 2)] + J. 

Consider next the matrix element Mi . We start with the following expression: 

1 „ 1 



Ml 



Gp 



dkk 



C- 



Hn-E k + Hi 



E 



H,p] 



(54) 



(55) 



Simple power counting shows that we can safely take at the origin, '5^1, and also replace the Green function G 
by its high-momentum asymptotics, G{p) — > — 2/^^. We then obtain symmetric and uniform expression 
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If we now rescale both p and p' as p ^ \j2k + 1 p, the integration over k factorizes and we obtain: 

Ml = -— [2/2(1, 1, 1) - 2/1(1, 1, 1) - /o(l, \f\ . (57) 

Finally, using explicit expressions for the integrals from Appendix, we arrive at our final result for the ultrasoft 
correction to the decay rate: 

161na2 /4 „, „ 20 \, 
Aus = 2 + ( — + 81n2+y-81nmjlna. (58) 



VIII. CONCLUSIONS 



The sum of all the contributions from Eqs. ( pl] , p"3|jl^j3^ , ^ , p8| ) gives the final result for the Oip? In a) corrections to 
Ps decay rate: 



.(0) 



p 

Ar, 



(0) ' TT 



3 r 



In^ a + Ina j --^ + 10 In 2 - 2Ap 
2 I 90 



— --ln^a + 7.9191na 

TT 2 



3, 2 
--In a 



, 229 Ao 
lna^-- + 81n2+- 



a 

TT 



In^ a. — 5.517 In a 

2 



(59) 
(60) 



Numerically, these corrections cause a negligible change in the theoretical prediction for p-Ps and o-Ps lifetimes at 
the current level of precision. It is interesting to note, however, that the magnitude of the leading 0(c? \\? a) and 
the subleading Oio? In a) corrections is comparable; in case of o-Ps they almost cancel each other. 

Our results Eqs.( [59|j60| ) are in agreement with two recent calculations of 0(p?\via) corrections In Ref. 

the result for 0(0;'^ In ct) correction to o-Ps decay rate has been obtained numerically, where as in Ref. [|j analytical 
methods similar to ours have been employed. We believe that the achieved agreement between three independent 
calculations ensures that the results, Eqs.( |59| , |60| ), are correct. 

As we mentioned, the Oict" In a) correction to Ps decay rates at present is not very interesting phenomenologically. 
A more important question, which we think we fully addressed in this paper, is how the logarithms of the fine structure 
constant can be efficiently extracted in the bound state calculation when the dimensional regularization is used to 
regulate the nonrelativistic dynamics. It is true that dimensional regularization offers many technical advantages in 
the calculation. This does not go without a price, however, since one has to be extremely careful in defining basic 
objects of the nonrelativistic theory, e.g. the wave functions and energies. If this is not done, one is left guessing 
whether or not is the calculation correct. 

Our key observation, which we think cures such problems and makes our calculation unambiguous, is the fact 
that the matrix elements in d dimensions are the uniform functions of the fine structure constant, and that the 
corresponding power of a can be determined by expressing the matrix elements in "d-dimensional" atomic units. We 
think that these arguments have not been spelled out before in the literature on one hand, and that they are necessary 
to make a convincing case, on the other. 

Finally, let us note that the technique discussed in this paper can obviously be used in other bound state QED 
problems, as well as for the heavy quarkonium states in QCD. 
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APPENDIX 



Definition of the integrals that were used in the derivation 

b 



V{a + h + c- d)T{a + b- d/2)r{b + c- d/2)T{d/2 - b) 
V{a)V{c)V{d/2)T{a + 2b + c - d) ' 



(62) 



/2(a, 6, c): 



d'^p' d'^p 
(27r)2rf 
T{a + b + c- 



= (47r 



d)r{a + b- d/2)V{d/2 - a)V{d/2 - b) 



r(a)r(6)r(c)r(d/2) 



(63) 
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